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Massive black hole binaries are the primary source of gravitational waves (GW) for the future 
eLISA observatory. The detection and parameter estimation of these sources to high redshift 
would provide invaluable information on the formation mechanisms of seed black holes, and 
on the evolution of massive black holes and their host galaxies through cosmic time. The 
Fisher information matrix has been the standard tool for GW parameter estimation in the last 
two decades. However, recent studies have questioned the validity of using the Fisher matrix 
approach. For example, the Fisher matrix approach sometimes predicts errors of > 100% in the 
estimation of parameters such as the luminosity distance and sky position. With advances in 
computing power, Bayesian inference is beginning to replace the Fisher matrix approximation 
in parameter estimation studies. In this work, we conduct a Bayesian inference analysis for 
120 sources situated at redshifts of between 0.1 < 2 < 13.2, and compare the results with 
those from a Fisher matrix analysis. The Fisher matrix results suggest that for this particular 
selection of sources, eLISA would be unable to localize sources at redshifts of 2 < 6. In contrast, 

Bayesian inference provides finite error estimations for all sources in the study, and shows 
that we can establish minimum closest distances for all sources. The study further predicts 
that we should be capable with eLISA, out to a redshift of at least 2 < 13, of predicting a 
maximum error in the chirp mass of < 1%, the reduced mass of < 20%, the time to coalescence of 2 
hours, and to a redshift of z ~ 5, the inclination of the source with a maximum error of ~ 60 degrees. 


PACS numbers: 04.30.Tv,95.85.Sz,98.90.Es,98.62.Js 


I. INTRODUCTION 

Through observations of nearby galaxies, it is appears 
that a supermassive black hole (SMBH) lies at the center 
of most, if not every, galaxy. Evidence further suggests 
that there is an intricate connection between the mass 
and evolution of the SMBH and the evolution of the host 
galaxy [1-3]. One of the major unsolved problems in 
astrophysics is how do the SMBHs form, and how do 
they affect galactic evolution over cosmic time? There 
are two main theories as to how seed black holes may have 
formed in the distant universe. In the first scenario, the 
first short-lived, low metallicity, stars produced low mass 
remnants at redshifts of 2 ~ 20. These so-called Pop III 
stars would have produced black hole seeds with masses 
of a few tens to a few thousand solar masses [4-11]. In 
the second scenario, the direct collapse of the centers of 
protogalactic disks would have produced black holes with 
masses of 10 5 — 10 6 M 0 at redshifts of 2 < 12 [12-15]. As 
of yet, there is no observational evidence to support one 
model over the other However, we do know that from 
observations of high redshift quasars, that SMBHs with 
masses of 10 9 M e existed at 2 ~ 7 [16]. 

It is also now accepted that whatever the mechanism 
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of seed formation, the most likely explanation for the 
masses of SMBHs that we observe in the local universe 
is a model of hierarchical structure formation through 
mergers [17-22]. While observational evidence strongly 
supports galactic mergers, the fate of the central SMBHs 
is unclear. It is believed that after the galactic merger has 
taken place, the two SMBHs sink to the bottom of the 
common potential well due to dynamical friction with 
the galactic environment [23, 24]. The two black holes 
then become gravitationally bound and form a binary 
(SMBHB) when the separation is on the order of < 10 
parsecs [25]. The subsequent evolution to the point where 
gravitational waves (GW) become the dominant mecha¬ 
nism of energy loss in the binary is still unknown, but 
is probably dependent on the role of mechanisms such 
as gas dynamics, the spins of the black holes, the stellar 
field near the galactic center and the possible existence 
of circum-black hole or circum-binary accretion disks. 

To compound the uncertainty in the final evolution 
of the SMBHB, there is little observational evidence for 
the existence of SMBHBs (see [26] for an interesting re¬ 
view). Among the current strongest candidates are: the 
radio galaxy 0402+379 [27] which has an estimated total 
mass of 1.5 x 10 8 Mq. The separation between the two 
SMBHs is estimated to be 7.3 pc, giving an orbital pe¬ 
riod of 150,000 years, and the distance to the source is 
2 = 0.055. OJ 287 is purported to contain two SMBHs 
with individual masses of (1.8 x 10 9 ,10 8 ) M 0 , and is at a 
redshift of 2 = 0.306 [28]. It is believed that the smaller 
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black hole orbits the larger object with an orbital pe¬ 
riod of 11-12 years. Finally, a potential sub-milliparsec 
binary was discovered in the galaxy SDSS J1201+30, at 
a distance of z = 0.146 [29]. Unfortunately, all of these 
candidates are at too low a redshift to usefully say any¬ 
thing about the high -2 behaviour of SMBHBs. Further¬ 
more, the errors in the mass estimates are on the scale of 
orders of magnitude, making it difficult to say anything 
useful on the black hole mass function in the universe. 

To try and answer some of these questions, ESA has 
recently chosen the theme of the “Gravitational Wave 
Universe” for the Cosmic Vision L3 mission selection. 
Within this program, a GW observatory called eLISA has 
been proposed [30, 31]. As the observatory will function 
in the frequency band of 10 -5 < f /Hz < 1, one of the 
main source targets will be the merger of SMBHBs to 
high redshift. We expect that with eLISA, it should be 
possible to measure the system parameters with sufficient 
accuracy that one could investigate the models of BH 
seed formation [32, 33]. 

In the last decade there has been a lot of progress made 
in the development of advanced algorithms that search 
the parameter space for SMBHBs. In general, because 
the signal-to-noise ratio (SNR) for these sources is so 
high, the actual detection is relatively simple. It is the 
estimation of the system parameters however that is the 
most time consuming. Until recently, the Fisher Infor¬ 
mation matrix (FIM) was the standard parameter esti¬ 
mation method in GW astronomy. The reason for this 
was two-fold : firstly, it was known from information the¬ 
ory that the inverse of the FIM provided the variance in 
the estimation of parameters. The fact that this calcu¬ 
lation could be done quickly (i.e. seconds to minutes), 
and with a relatively small investment in coding, made 
a FIM a popular choice. Secondly, in the late 1980s, 
the emerging field of information geometry demonstrated 
that the FIM could be associated with a metric tensor on 
the parameter manifold, and to the local curvature of the 
log-likelihood. At a time when template banks were the 
only conceivable option to search for GWs, the FIM was 
implemented in the proper spacing of templates in the 
bank [34-37]. The early success of the FIM ensured that 
it stayed as the main parameter estimation tool within 
the GW community. 

In the last decade however, and mostly due to the leap 
in computational power, there has been a paradigm shift 
within the community when it comes to parameter esti¬ 
mation. It was demonstrated very early on that template 
banks would be computationally prohibitive for a space- 
based observatory [38-40]. This pushed the space-based 
GW community to move to more dynamic Bayesian 
based algorithms. For the question of parameter es¬ 
timation, different variants of a Markov Chain Monte 
Carlo (MCMC) became very popular [40-52]. While the 
method required a much higher investment in implemen¬ 
tation of the algorithm, and while it took much longer to 
run (hours to days), the benefit was that one obtained the 
marginalized posterior distributions for each of the pa¬ 


rameters and allowed for the easy incorporation of prior 
information. If the algorithm had converged, then one 
had all the information that one could ask for regarding 
the problem at hand. 

However, what began as an effort to implement an al¬ 
ternative method that would confirm and solidify the re¬ 
sults of the FIM, we very quickly began to observe sit¬ 
uations where the results from the two methods did not 
agree. On further investigation, we began to see many 
problems with the FIM calculation (e.g. 100% errors in 
the estimation of the luminosity distance, or errors in 
the sky position larger than the entire sky). For the first 
time, the community began to question the validity of 
using the FIM (see, for example, Ref [53] for an in-depth 
analysis). As time has progressed, confidence in the reli¬ 
ability of Bayesian inference methods have grown, while 
the FIM is now viewed as an approximation method to 
be used with caution. 

All of this has led to an unfortunate, and unpalatable 
situation, where we are now using advanced Bayesian 
inference methods to demonstrate parameter estima¬ 
tion capabilities for individual sources, of many different 
source types, but revert to the FIM when large scale pa¬ 
rameter estimation studies need to be conducted, even 
though we know a-priori that many of the results will be 
incorrect. Unfortunately, these types of studies are usu¬ 
ally associated with science impact investigations for mis¬ 
sion formulation, upon which policy decisions are made 
(see, for example [31, 54]). As, to our knowledge, this 
is the largest ever Bayesian inference study undertaken 
for SMBHBs, the goal of this work is to take the first 
step towards large-scale Bayesian inference studies, and 
to caution against drawing conclusion from a FIM study. 


A. Outline of the paper 

The paper is structured as follows. In Sec II we define 
the eLISA mission concept, and outline the response of 
the observatory to an impinging GW. We also define the 
higher harmonic inspiral waveform that will be used for 
the study. Sec. Ill introduces the Fisher information ma¬ 
trix and discusses potential failings of the method. This 
is followed by a description of the MCMC algorithm im¬ 
plemented in this work. In Sec. IV we discuss the gen¬ 
eration of the sources that were studied, as well as the 
parameter priors used for the Bayesian inference. Sec. V 
presents the main results comparing the Fisher matrix 
and Bayesian inference methods. 

II. THE eLISA DETECTOR RESPONSE 

The European Space Agency (ESA) has recently cho¬ 
sen the theme of the “Gravitational Wave Universe” for 
the L3 mission within the Cosmic Vision framework. For 
this theme, a mission concept called eLISA has been pro¬ 
posed. This mission is a restructured version of the LISA 
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good approximation for the detector response up to a fre¬ 
quency of 0.2 Hz. It was demonstrated that the inclusion 
of higher harmonic corrections into the comparable mass 
waveform dramatically improves parameter estimation 
and reduces correlations between parameters [56-60]. It 
is clear that with reduced sensitivity at low frequencies, 
the higher harmonic corrections will be indispensible for 
detecting and resolving the most massive binaries. With 
this in mind, the post-Newtonian (PN) inspiral waveform 
polarizations, including higher harmonic corrections up 
to 2-PN order are defined by [61] 


H-,x — 


( 3 ) 


FIG. 1: Instrumental noise model for a 10 9 m arm eLISA con¬ 
figuration. 


mission, composed again of a constellation of three satel¬ 
lites, following ballistic heliocentric orbits. The main dif¬ 
ferences between eLISA and LISA are (i) a shortening of 
the arm-length between spacecraft from five to one mil¬ 
lions kilometers, and (ii) four, rather than six, laser links. 
The idea behind the mission is to launch the three satel¬ 
lites into an orbit roughly 20° behind the earth, inclined 
at 60° to the ecliptic. The constellation is composed of 
one “mother” spacecraft at the center of the constella¬ 
tion, with two outlier “daughter” spacecraft. The con¬ 
stellation retains a 60° angle between spacecraft, with a 
nominal mission lifetime of two years. In Figure 1 we 
present the eLISA instrumental noise curve. Compared 
to the LISA noise curve, we are worse off in sensitivity 
by a factor of 5 at lower frequencies, and the “bucket” 
of the noise curve has moved upwards and two the right. 
This makes the constellation more sensitive to lower mass 
binaries. At the highest frequencies, we gain a little sen¬ 
sitivity due to a reduction of photon shot noise, coming 
from the shorter armlengths. 

In the low frequency approximation [55], where the 
GW wavelength is greater than the arm length of the 
detector, the strain of the gravitational wave (GW) in 
the eLISA detector with both polarizations is given by 
the linear combination 

h(t) = h + m)F + +h><m)F x • (i) 

Here we define the phase shifted time parameter £(t) as 

£(f) = t — I?® sin d cos (a(t) — </>), (2) 

where R® = 1 AU /c is the radial distance to the detector 
guiding center, c is the speed of light, ( 6 , </>) are the po¬ 
sition angles of the source in the sky, a(t) = 27 xf m t + k, 
f m = 1 /year is the LISA modulation frequency and k 
gives the initial ecliptic longitude of the guiding center. 
Due to the shorter arm length of eLISA, the LFA is a 


Here m = nil + m 2 is the total mass of the binary, 
V = mi m 2 /to 2 is the reduced mass ratio and Dl is the 
luminosity distance of the source. We can define a re¬ 
lation between redshift, z, and luminosity distance, Dl, 
within a ACDM model by 

D — (1 ) ° f Z ^ Z ' 

F 0 Jo + z ; ) 4 + Hm(1 + z') 3 -t- Ha 

( 4 ) 

where we use the concurrent PLANCK values of fin = 
4.9 x 10 -5 , JIM = 0.3086 and = 0.6914 and a Hub¬ 
ble’s constant of Hq = 67 km/s/Mpc [62], The invariant 

PN velocity parameter is defined by, x = (Gmuj/c 3 ) 2 ^, 
where w = d^orb/dt is the 2 PN order circular orbital 
frequency and 4> orf) = 1/2 <&gw — <P% rb ~ <fiorb(.t) is the 
orbital phase defined by 
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where 0 is related to the time to coalescence of the wave, 
t c , by 

m = Sk {tc ~ th (6) 

and ip° rb is the orbital phase of the wave at coalescence. 
For the rest of the paper, we will work with GW phases. 
In Eqn (3), the functions contain the PN correc¬ 

tions to the amplitude and the phase harmonics. All 
of the half-integer, and some parts of the higher integer 
1 T+ x terms contain the factor Sm = mi — m 2 which 
supress the odd phase harmonics and parts of the even 
phase harmonics in the equal mass case. Furthermore, 
we see the extra frequency harmonics arising due to the 
higher order x n terms, and the fact that the different 
H+ x terms contain different multiples of n& or t,. 

The functions F +,x are the beam pattern functions of 
the detector, given in the low frequency approximation 
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by 


Fk{t-,ip,9,(/>) = 

+ 


\ [cos(2i p)D + (t; ip, 6, 0, X k ) 
sin(2 il>)D x (t; if), 6 , (j>, A fc )] , 
i [sm(2ip)D + (t-,il),0,<j), X k ) 
cos(2ijj)D x (t; ip, 6, <j>, A*,)] , 


(7) 

( 8 ) 


where ij) is the polarization angle of the wave and we 
take A = 0 for the single channel eLISA configuration. 
Explicit expression for the detector pattern functions 
D +,x (t) can be found in [63]. The above equations for the 
frequency and phase can sometimes break down before 
we reach the last stable circular orbit (LSO) at R = 6 M. 
To compensate for this, we terminate the waveforms at 
R = 7 M in the cases where the coalescence time is 
shorter that the duration of the signal, t c < T 0 t s . For 
low mass sources where the LSO frequency is possibly 
close to the limit or outside the band of the detector, 
we terminate the waveform when the frequency of the 
highest harmonic reaches 0.2 Hz. 


III. STATISTICAL METHODS FOR 
ESTIMATING PARAMETER ERRORS. 

A. The Fisher information matrix. 


One of the main tools used in the GW community for 
the estimation of parameter errors is the Fisher infor¬ 
mation matrix (FIM). In the high SNR limit, assuming 
that the waveform is a linear function of the waveform 
parameters, and the error distributions in the parameter 
errors are Gaussian, the inverse of the FIM is said to give 
the variance-covariance matrix. The square root of the 
diagonal elements of the inverse FIM thus give a 1-a es¬ 
timation of the error in our parameter estimation. Given 
a template h(t\ A K ), the FIM is defined by 
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dh dh 
5W 5A 7 


= -E 


d 2 In C 
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where h = h(f; A K ) is the waveform template in the 
Fourier domain and the angular brackets denote the 
Fourier domain noise-weighted inner product 


(h\s) = 2 



df 

Sn(f) 




( 10 ) 


where a tilde denotes a Fourier transform and an asterisk 
denotes a complex conjugate. The quantity S n (f) is the 
one-sided noise spectral density of the detector, which 
is normally a combination of instrumental and galactic 
confusion noise, but in this article represents the instru¬ 
mental noise model for eLISA presented above. Finally, 
the FIM can also be interpreted as the second derivative 


of the log-likelihood, ln£, where we define the likelihood 
given a detector output s(t) and a GW template h(t) as 

£(A«) = exp (-i {s - h(X K ) \s - h( A*) >) . (11) 

For almost two decades, the FIM has been the primary 
tool within the GW community when it comes to param¬ 
eter estimation. It has mostly gained popularity from 
studies on the construction of template banks as it was 
shown that there is a close relationship between the FIM 
and the metric tensor on the parameter search space [35— 
37]. This dual interpretation of the FIM as both a sta¬ 
tistical and geometric tool have proved to be very useful. 
A further advantage of the FIM is the speed of calcula¬ 
tion. In general, one could calculate the FIM in seconds 
to minutes, depending on the mass of the system and 
the duration of the signal. This ability to obtain a quick 
error estimate for the system parameters concretized the 
usefulness of the FIM for GW astronomy. 

However, the FIM has not been without its difficul¬ 
ties. The first, and probably the most influential, is the 
interpretation of the FIM results. This depends heavily 
on whether the problem in hand is attacked from a fre¬ 
quent ist or Bayesian approach. Ref [53] has very nicely 
laid out the main ways of interpreting the FIM, which 
also clearly displays the confusion arising in its interpre¬ 
tation : (i) is a lower bound (Cramer-Rao) for the 
error covariance of any unbiased estimator of the true 
parameter values. It is therefore a frequentist error for 
an experiment characterized by a true signal and a par¬ 
ticular noise realization, (ii) is the frequentist error 
in the maximum likelihood estimator for the true param¬ 
eter, assuming Gaussian noise and a high signal to noise 
ratio (SNR), (iii) gives the covariance of the pos¬ 
terior probability distribution for the true parameters, 
inferred from a single experiment, assuming the presence 
of a true signal, Gaussian noise and a high SNR. In this 
respect, the FIM provides a level of Bayesian uncertainty 
rather than an error. 

As most GW transient sources fall into the final cat¬ 
egory and require a Bayesian interpretation, one would 
normally use interpretation (iii) for the FIM. However, 
this interpretation also has a number of other problems. 
Firstly, as can be seen from Eq. 9, the FIM is a pure 
template calculation which assumes a certain model for 
S n (f). The actual detector data is not used for the cal¬ 
culation. As we can also see from the right hand side of 
the same equation, the FIM can be further interpreted as 
the expectation value for the local curvature of the log- 
likelihood. In this respect, the FIM has no knowledge 
of other possible modes of the solution or on the general 
structure of the likelihood surface. Finally, even if the 
SNR is large, V KU can be either singular, or at least ill 
conditioned (which happens quite frequently in GW as¬ 
tronomy), which results in problems as r K „ has in general 
to be inverted numerically. The only possible solution in 
this case is then a more robust Monte Carlo study. 

Another concern with the FIM comes from the very 
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statistics-geometry duality that established its place in 
the GW community in the first place. Geometrically, the 
FIM can be interpreted as being related to both the lo¬ 
cal metric tensor on the parameter space, and as a local 
estimate of the curvature of the log likelihood function. 
Crucially however, the FIM is not a coordinate in¬ 
variant quantity. This immediately poses the problem 
of what is the best choice of coordinates for a particular 
problem. For example, it was demonstrated in Ref [64] 
that one could approximate the true likelihood for some 
sources using the FIM. It was further demonstrated that 
using a FIM with a u = 1/Dl parameterization provides 
a much closer approximation to the true likelihood than 
a parameterization in Dl itself. Furthermore, as stated 
earlier, one of the primary conditions for using the FIM 
in a statistical analysis is that the distributions in pa¬ 
rameter errors are Gaussian. In reality, data in physics, 
astrophysics, economics etc. is rarely normal. In this 
and other studies, we have observed a number of situa¬ 
tions where the posterior distribution in Dl, for example, 
follows a log-normal distribution. This suggests that for 
these sources, an error estimate from a FIM using a In Dl 
parameterization will be more exact than an estimation 
using Dl- However, this information is obtained via an a- 
posteriori comparison with another method, so no global 
statement can be made on the a-priori choice of parame¬ 
terization for a problem. We should also again highlight 
the fact that in GW astronomy V KU is almost always ill- 
conditioned. In this case, working with the bare physical 
parameters returns a condition number for T KV that effec¬ 
tively prevents a numerical inversion of the matrix. How¬ 
ever, working in logarithms of the very small/large val¬ 
ued parameters sufficiently reduces the condition number 
such that in most, but not all cases, we can numerically 
invert T KV . 

To further complicate matters, there is also an impor¬ 
tant interplay between the geometry-statistics duality, 
and the physical problem at hand which is often over¬ 
looked. If we take the GW problem that we are study¬ 
ing as an example, we can formally state that the source 
search and resolution takes place on a n-dimensional Rie- 
mann manifold, with well defined norm, scalar product 
and metric tensor. We could initially define all points on 
this manifold as A K £ R n . One could further state that 
the physical parameter set is defined by the following 
subsets : {mi, m2} £ [10 2 ,10 s ] M@, z £ [0, 20], t c £ [0, 2] 
years, {9,i,if} £ [0, 7 r] and {</>, ip c } £ [0, 2n]. To high¬ 
light the problem, let us start with the R 2 sub-space 
defined by the mass parameters. As well as the coordi¬ 
nates (mi, 777-2), one could also construct a sub-manifold 
by combining any of the two parameters {m, 77 , q , M c , 77 }, 
where q = mi/m 2 and all other quantities have been 
previously defined. It is well known that the coordinates 
(m\, m2) are not a good choice as the parameter space is 
degenerate, and any map of the form (m, if) —> (mi, m 2 ), 
for example, is not a diffeomorphism due to the map not 
being bijective. This immediately constrains the choice 
of “good” mass parameters. So, for arguments sake, let’s 


say that we choose to work with the coordinates (q,if), 
and let us further complicate matters by assuming that 
we are investigating an almost equal mass binary. One 
could imagine that our robust detection technique pro¬ 
vides solutions of q = 1.1 and 77 = 0.2494 for this partic¬ 
ular source. Our subsequent FIM calculation then pro¬ 
vides a ±1 a error in each of these parameters. However, 
for binary systems q m in = 1 and r] max = 0.25, so what 
exactly does it mean if the FIM predicts error bound val¬ 
ues lower/higher than these values? In this case, the FIM 
is blameless. It assumes a Gaussian distribution of errors 
and is unaware of any boundaries imposed by the physi¬ 
cal problem at hand. Geometrically, it assumes that the 
(q, rf) sub-manifold is infinite and continuous. So while 
there are no geometrical/statistical problems with this 
choice of coordinates, astrophysics places a finite bound¬ 
ary on the manifold. It would therefore seem that for 
compact binary systems, the best choice of mass param¬ 
eters is any combination of {m,M c , n}, as long as the 
mapping to (mi, m 2 ) space is diffeomorphic. 

A final comment involves the structure of the manifold 
and the physical boundaries on the angular parameters 
{6, (f), if, ipc}. We stated earlier that our problem was one 
which could be formulated in R n . More correctly, due to 
the sky position coordinates, our manifold actually has a 
S 2 x R n ~ 2 structure, i.e. the merging of a 2-sphere with 
a R m subspace, where m = n — 2. Experience seems to 
suggest that the FIM always “assumes” we are working 
on an infinite R n manifold with no boundaries. The fact 
that the 2-sphere S 2 is finite may have a knock-on ef¬ 
fect on the FIM calculation, and may be why the FIM 
sometimes returns error estimates that are larger than 
the area of the 2-sphere. Furthermore, if we examine the 
response function, it should be clear that h(t) is invariant 
under a translation of the form (ip — > if ± n, ip c —> ± 27 t ), 
justifying the physical boundaries of 7 r and 27 t defined 
above. However, at times the FIM can predict 1 <t error 
estimates in these parameters that are many multiples of 
7 T in size. Common practice would then be to say that 
the error on these parameters were the size of the phys¬ 
ical boundary. However, it was shown in Ref [40] that 
doing so causes an underestimation of some of the other 
parameters. It was found that, for example, the (M c , ip c ) 
sub-space is a cylinder embedded in the larger param¬ 
eter space. On this sub-manifold, the likelihood peak 
was wrapped multiple times around the cylinder, with a 
length that was in accordance with the FIM prediction. 

All of this demonstrates that if one is truly intent on 
using the FIM, one must (a) make sure that the choice 
of coordinates is optimal, (b) ensure that any mappings 
between the FIM coordinates and true physical parame¬ 
ters are diffeomorphic, (c) make sure that the geometry- 
statistics duality is fully compatible with the astrophysi- 
cal implications , (d) make sure that the condition num¬ 
ber is sufficiently small that the results of a numerical 
inversion are believable and (e) even if everything seems 
above board, and all constraints have been respected, be 
very careful with the interpretation as it may still be in- 
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correct. 


B. Markov Chain Monte Carlo 


The goal of any Bayesian analysis is to test a model 
hypothesis for a particular problem. For GW astronomy, 
given a data set s(t), a waveform model based on a set 
of parameters A K and a noise model for S n (f), one can 
construct the posterior probability distribution p(A K |s) 
via Bayes’ theorem 


p(A K |s) 


^(A*>(s|A K ) 

p(s) 


( 12 ) 


The prior probability 7 t(A k ) reflects of assumed knowl¬ 
edge of the experiment before we evaluate the data. 
p(s|A K ) = £(A K ) is the likelihood function defined by 
Eq. 11, and p(s) is the marginalized likelihood or “evi¬ 
dence” given by 


P(s) = J d\ K n(X K )p(s\X K ). (13) 


A common method for carrying out Bayesian analysis is 
to use a Metropolis-Hastings variant of the MCMC al¬ 
gorithm [65, 66]. This quite simple algorithm works as 
follows : starting with a signal s(t) and some initial tem¬ 
plate h{t\ A K ), we choose a starting point x(A K ) randomly 
within a region of the parameter space bounded by the 
prior probabilities 7r(A K ), which we assume to contain the 
true solution. We then draw from a proposal distribution 
and propose a jump to another point in the space x'{X v ). 
To compare the performance of both points, we evaluate 
the Metropolis-Hastings ratio 

H _ TT(x')p(s\x')q{x\x') ^ 

n(x)p(s\x)q(x'\x) 

Here q{x\x') is a transition kernel that defines the jump 
proposal distribution, and all other quantities are previ¬ 
ously defined. This jump is then accepted with probabil¬ 
ity a = otherwise the chain stays at x(X K ). 

In order to improve the overall acceptance rate, the most 
efficient proposal distribution to use for jumps in the pa¬ 
rameter space is, in general, a multi-variate Gaussian dis¬ 
tribution. As we assume, in most cases, that the FIM will 
reliably describe the likelihood surface for the problem at 
hand, a good starting point is to thus use the FIM. The 
multi-variate jumps use a product of normal distribu¬ 
tions in each eigendirection of r K „. The standard devi¬ 
ation in each eigendirection is given by a K = 1/y/DE K , 
where D is the dimensionality of the search space (in 
this case D = 9), E K is the corresponding eigenvalue of 
V Kl/ and the factor of 1 /y/D ensures an average jump of 
~ lcr. In this case, the jump in each parameter is given 
by <5A K = Af( 0,1 )a K . This type of MCMC algorithm is 
commonly referred to as a Hessian MCMC and is known 
to have the highest acceptance rate for this particular 
family of random walk algorithms. 


One problem with the MCMC approach is that, while 
it provides the full marginalized posterior density func¬ 
tion (pdf) for each of the system parameters, it is essen¬ 
tially a black-box algorithm. The MCMC provides reli¬ 
able answers provided the chain is run for long enough, 
and is efficiently sampling the parameter space. The 
idea of what constitutes “long enough” is both appli¬ 
cation, and in our situation, source dependent. Thus 
MCMC algorithms have to be tailored to the problem in 
hand. While there are some defined convergence tests 
for MCMC algorithms, they are usually tenuous at best. 
The best convergence tests are still “by-eye” tests, such 
as visually examining the convergence of the means and 
standard deviations of the chains. The main issue then 
becomes ensuring long enough chains to have confidence 
in our predictions. This means that MCMC algorithms 
take time to map out the pdfs, making a quick estimation 
out of the question. 

As the convergence of the MCMC can be slow, espe¬ 
cially for low mass, high z and/or low SNR systems, we 
use a mixture model combining a Hessian and a Differ¬ 
ential Evolution (DE) MCMC to estimate our parame¬ 
ters (hereafter DEMC) [67, 68 ]. Differential Evolution 
is a subset of a genetic algorithm that normally works 
by combining information from multiple simultaneous 
chains and using the information to construct jump pro¬ 
posals. However, the standard implementation of DE 
can be computationally prohibitive, especially for long 
waveform lengths. To overcome this problem, we use a 
modified version of the DE algorithm. Once the burn-in 
part of the chain is over, we build a trimmed history of 
the chain, i.e. we keep every 10th chain point. For each 
parameter, this then means that we have the full chain 
array Xi of instantaneous length n, and a trimmed history 
array Xi with instantaneous length of n = nj 10. We then 
construct our proposed point in parameter space using 

x i+ i = Xi +-y(xj - x k ) , (15) 

with j , k £ U[ 1, n] and i ^ j ^ k. The parameter 7 has 
an optimal value of 7 = 2.38/V2 D, where again D is the 
dimensionality of the parameter space. In this implemen¬ 
tation of the algorithm, we use a (2/3, 1/3) split between 
DE and Hessian MCMC moves. Furthermore, we know 
that due to the shape of the beam pattern functions, the 
response of an observatory to an impinging GW is natu¬ 
rally bimodal. To account for this, we set 7 = 1 for 10% 
of the DE jumps. This helps moves the chain between 
modes of the solution and encourages greater exploration. 

Finally, whether a move is proposed by the Hessian 
or DE part of the algorithm, the choice to move or 
stay requires a calculation of the likelihood. To ac¬ 
celerate this calculation and to allow us to investigate 
the lower mass and/or high redshift systems, we use 
the composite integral method described in Ref [69] to 
speed up the computational run-time of the algorithm. 
To aid practicality for the problem of non-spinning 
SMBHBs, the parameter set we chose to work with is 
A K = (ln(M c ), ln(/z), ln(t c ), cos 9, <fi, ln(Z)/,), cos t, <p c , ip} 
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FIG. 2: eLISA redshift horizon distance (excluding merger 
and ringdown) as a function of redshifted total mass for mas¬ 
sive black hole binaries assuming a threshold of pth = 10 on 
the inspiral-only SNR. 


where M c = mrf / 5 is the chirp mass, p, = mrj is the 
reduced mass and all other parameters have been previ¬ 
ously defined. 


IV. ASTROPHYSICAL SOURCE SELECTION 

In order to select the specific test-sources for this 
study, we first generated a population of SMBHBs from a 
Monte Carlo simulation based on an astrophysical source 
model [19]. To ensure that the comparable mass PN 
waveform was still valid for the types of sources com¬ 
ing from the population model, we only accepted sources 
with mass ratios of q £ [1,30]. Furthermore, to allow 
us to run long-enough chains in order to have confidence 
in our Bayesian inference, we restricted all waveforms to 
an array length of n < 2 21 elements. This ensured that 
our Markov chains would run in an acceptable amount 
of time. While the mass, mass ratio and redshifts came 
from the distribution, other constraints were imposed by 
hand : the time to coalescence for each source was ran¬ 
domly chosen between t c £ [0.3,1] yr and all angular 
values were chosen over their physical ranges. To choose 
a detection threshold, it was recently shown [70] that a 
null test (i.e. the detector output consists of instrumen¬ 
tal noise only : s(t) = n{t )) returns a detection SNR 
threshold of pth = 10 for the eLISA configuration, where 
we define the SNR as 


V( h \ h ) 


(16) 


This value was again imposed for this study. In Fig 2 we 
plot the detection horizon z for the sources that returned 
P > pth as a function of total redshifted mass. We should 
mention here that the detection horizon is slightly larger 


to the one presented in Ref [70] as a more extensive Monte 
Carlo was run for this study. 

As the goal of this work is the comparison between dif¬ 
ferent parameter estimation methods, and not the astro- 
physical significance of source detection and resolution, 
we then trimmed the the generated population to the set 
of study sources as follows : our sources covered a red- 
shift range of 0.1 < z < 13.25. We then divided this 
range into redshift bins of A z = 0.1, and from each bin 
drew a single source, providing us with 120 sources. In 
Fig. 3 we plot the total redshifted masses m(z) of the 
selected sources as a function of redshift in the left hand 
panel, and their optimal SNRs as a function of redshift in 
the right hand panel (i.e. the SNRs obtained by using the 
true parameter values). We can see that from redshifts 
of z > 6 onwards, we very quickly begin to approach the 
detection threshold. While two of our sources fell just 
below the detection threshold, we decided to study them 
anyway to investigate if we could detect and resolve the 
sources. In these cases, our chains actually found solu¬ 
tions that differed from the true solution (as would be 
expected due to the presence of noise), but with SNRs 
that were superior to the detection threshold. We should 
point out here that in many of the redshift bins, brighter 
sources existed. As a consequence, one should be careful 
in interpreting the results here as an absolute limit of ca¬ 
pability for the eLISA detector. In fact, due to the pres¬ 
ence of brighter sources in the redshift bins, for many of 
the sources presented here, the results can be interpreted 
as a “worst-case” scenario. 

Finally, for each source we ran a 10 6 iteration DEMC, 
with a “burn-in” phase of 2 x 10 4 iterations. During the 
burn-in phase, we used a combination of thermostated 
and simulated annealing to accelerate the chain mixing 
and convergence to the global solution [40]. With each 
chain we checked the instantaneous means and standard 
deviations to ensure they had converged to constant val¬ 
ues. One thing that we observed, was once again how 
important it is to have long chains to properly attain 
convergence, especially for low mass and/or high redshift 
sources. 

One further note is that for z > 3, the posterior densi¬ 
ties for the vast majority of the sources become at least 
bi-modal, if not multi-modal. It has been well docu¬ 
mented that the detector response h(t) is invariant un¬ 
der certain symmetries, either in the polarizations of 
the waveform, or in the beam pattern functions. For 
example, a change such as (9, <j>) — > (n — 9, <j> ± 7r) or 
(0, <f>,L,ip) —I (tt — 9, <f> ± 7T, 7r — i, 7r — ip) leaves h(t) in¬ 
variant. If the signal is strong enough, the Markov chains 
normally stay on one mode, as it is to unlikely within the 
run-time of the algorithm to walk across the parameter 
space to the other mode(s). In general, if many chains 
are run simultaneously, we would expect a 50-50 split be¬ 
tween chains ending up on the main mode, and say, the 
mode corresponding to the antipodal solution in the sky. 
It was clear however, from initial runs, that for systems 
where the sky error is large, the likelihood surface is flat 
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FIG. 3: The left hand cell plots the total redshifted mass m(z ) of the 120 selected systems as a function of distance. The right 
hand cell displays the optimal SNR of each source, also as a function of distance. 


enough that the chain does not have to move very far in 
parameter space before it can change modes. In these 
cases, we ran extended chains and evaluated results for 
the dominant mode solution. 

We chose to use flat priors in 
{In Dl, In M c , In /i, lnt c , cos l , cos 9, </>, ip, ip c }, where 
the boundaries were set at Dl E [7.7 x 10~ 4 ,300] Gpc, 
M c G [435,5.06 x 1O 7 ]M 0 , p, G [250,2.9 x 10 7 ]] M 0 , 
t c G [0.2,1.1] years, {cos i, cos 9, <p} G [—1,1] and 
were left completely open for {ip,ip c }. The prior in 
luminosity distance is chosen such that we assume that 
there are no SMBHBs closer than the M31 galaxy, and 
no further than a redshift of z ~ 25. The priors in 
{M c , p} are chosen so that the minimum total redshifted 
mass is m(z) = 1O 3 M 0 , while the upper limit of 
m(z) = 1.163 x 10 8 M 0 assures that the maximum 2PN 
harmonic last stable orbit frequency corresponds to a 
value of 5 x 10 -5 Hz. These values also correspond to 
a prior range in reduced mass of 1 < q < 100. The 
choice of priors that are flat in the logarithm of the 
dimensionful quantities {In Dl, In M c , In fi, In t c } says 
that these parameters are unknown a-priori to within 
an order of magnitude across a wide prior range. In 
practice a better approach would be to adopt an astro- 
physically motivated joint prior on {ln.Dj,,hiM c ,ln/x} 
that takes into account the number of mergers as a 
function of redshift and the mass evolution of the 
population. By assigning hyper-parameters to such a 
prior and considering the joint posterior distributions 
for all SMBHBs detected by eLISA, it will be possible to 
constrain the merger model, in much the same way that 
the distribution of white dwarf binaries in our galaxy 
can be contained by eLISA [64]. We leave such a study 
to future work, but note for now that a uniform prior 


in In Dl corresponds to a 1 /Dl prior on the luminosity 
distance Dl, while the FIM study implicitly assumes a 
uniform prior in 1 /Dl, which corresponds to a 1 /D\ 
prior on Dl- This should be contrasted to a more 
astrophysically motivated prior that grows as D\ for 
small Dl, before rolling off to some decaying function of 
Dl for large luminosity distances. Since the results will 
be dominated by the likelihood at small Dl and by the 
prior at large Dl, our choice of priors that decay in Dl 
are not unreasonable. 


V. RESULTS 

To compare the FIM and a full Bayesian inference, our 
plan of action is as follows : the first step is to test the 
hypothesis of Gaussianity of the parameter posterior den¬ 
sity distributions, and as a consequence, the validity of 
using the FIM for parameter estimation. Once we obtain 
the levels of Gaussianity, we calculate the FIM at the 
median of the distribution, X K , for each of the parame¬ 
ters { M c , p, t c , Dl, t}- As we are carrying out a Bayesian 
inference study, our goal is to estimate the 95% Bayesian 
credible intervals (BCI) for each of the parameters. For 
the FIM estimate, we express the 95% quantiles QI , as 
calculated at X K , and assuming a Gaussian distribution 
in the parameter errors, are simply given by 

QI 95% = X K ±lMa K , (17) 

where cr K is the standard deviation of the distribution 
obtained from a K = T KK ' . As shown earlier, Bayes’ 
theorem can be written in the form 

p (A K |s) oc 7T (A K ) C (A K ). 


(18) 
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FIG. 4: A plot of skewness versus excess kurtosis for each of the 120 sources in the study. The results display that the vast 
majority of the sources in the study have non-Gaussian posterior density distributions. 


Once we have the posterior distribution p (A K |s) from the 
chains, we can now calculate a credible interval C such 
that 


J p (A K |s) dX K = 1 — a, 


(19) 


where for a 95% BCI, a = 0.05. This allows us to make 
a degree-of-belief statement that the probability of the 
true parameter value lying within the credible interval is 
95%, i.e. 


Wheels) =0.95. (20) 

For the positional resolution of the source, we can define 
an error box in the sky according to [55] 

Afl = 27ry / £ ee £«’ - (E^) 2 , (21) 

where 

E ee = (A cos 9A cos 0 ), (22) 

E^ = (A0A 0), (23) 

£** = (Acos0A<p) , (24) 


and E KI/ = ^AA fc AA") are elements of the variance- 
covariance matrix, found by either inverting the FIM, 
or are calculated directly from the chains themselves. 


A. Testing the hypothesis of Gaussianity 


The standard assumptions for the validity of using the 
FIM are : we in the high SNR limit, and we assume that 
the errors in the estimation of the system parameters are 
described by a multivariate Gaussian probability distri¬ 
bution 


P( AA K ) 



-ir^AA^AA" 


where, once again, we define the FIM 


r 


K1S 


dh dh 


= -E 


' a 2 In £ ' 

d\ K d\" 


(25) 


(26) 


as the negative expectation value of the Hessian of the 
log-likelihood, and |F| = det(Y KV ). To test these hy¬ 
potheses, we calculated the skewness and kurtosis of the 
posterior density distributions p(X K \s). 

Using the chain points we can now investigate the devi¬ 
ation of each parameter distribution from a normal dis¬ 
tribution. As is standard, we begin with the fact that 
a random variable x is said to be normal if its density 
function follows the form 


f(x) 


V2^o 




(27) 


where p and a are the mean and standard deviation re¬ 
spectively. From this we define the skewness and excess 
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FIG. 5: In the left hand cell we represent the 95% quantile intervals for the luminosity distance based on a FIM calculation 
at the chain median. The right hand cell represents the 95% BCI for a Bayesian inference analysis. Not only do the FIM 
intervals not contain the true value (represented by the dashed line) for some sources, but for many sources at z > 6 suggest 
that the source is unresolvable with a > 100% error in the distance estimate. However, a Bayesian inference shows that the 
FIM estimation is incorrect and all intervals are not only finite, but importantly, contain the true value. 


kurtosis via the third and fourth standardised moments 


\fl i 
^2 


E(x — p) 3 E(x - p) 3 
[E{x - pi) 2 } 3 / 2 = V 3 ’ 

E(x~p) 4 0 _E(x-p) 3 

[E(x - p) 2 } 2 a 4 


(28) 

(29) 


where again E denotes the expectation value. We note 
here that we are using the excess kurtosis (i.e. with the 
factor of -3) as for a normal distribution, y/fii = @2 = 
0. For historical reasons, the skewness is defined in this 
manner even though we can have \//3i < 0. For a sample 
of n data points, i.e. xi,..,x n , we can write the the 
sample estimates of yffti and f3 2 in the form 


{M c ,p,t c , Dl,o}- The results are plotted in Fig. 4. We 
remind the reader that if the Gaussian assumption of er¬ 
ror distribution was true, and hence it was valid to use 
the FIM, both \fB\ and f3 2 would be clustered around 
zero. It is clear from the figure that the vast majority 
of the marginalized posteriors fail the Gaussianity test 
for all of the parameters. This strongly suggests that 
we should not believe the results of the FIM for these 
sources and should be encouraged to use other methods 
when carrying out a full parameter estimation for SMB- 
HBs. 


B. Luminosity Distance 


Vh = ^ (30) 

m 2 


where each moment m k is defined by 

m k = ^2 {xi - x) k /n (32) 

i 

and the sample mean x is given by 

x = ^ Xi/n. (33) 

For each of the 120 sources, we analysed the skewness 
and excess kurtosis of the marginalized posterior densi¬ 
ties for the most astrophysically interesting parameters 


As stated in the introduction, the main driving force 
behind this study are previous works where results based 
on the FIM suggest that A Dl/Dl > 1. For these sys¬ 
tems, the results imply that we would be unable to locate 
them anywhere in the Universe. Conceptually, however, 
this is clearly untrue. As an example, it may be that we 
detect a source at z = 10 with p = 11. A FIM analysis 
returns a value of A Dl/Dl > 1 for this source, giving a 
100% or greater error in the distance estimation. How¬ 
ever, if we were able to transport that source from z = 10 
to somewhere in the local group of galaxies, it would have 
an SNR in the millions. The fact that observed SNR is 
small clearly implies that this source is not that close to 
us, and in fact, cannot be closer than a certain minimum 
distance. 

In Fig. 5 we plot the 95% QI for the luminosity distance 
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FIG. 6: The left hand cell represents the sky error box in square degrees based on the FIM for the 120 selected systems as a 
function of distance, while the right hand cell plots the Bayesian inference error box for each source. The upper limit of each 
cell corresponds to the total area of the sky. We can see that for sources at z > 9, the FIM results suggest that the error in sky 
location can be equal to or greater than the size of the total sky. However, we can see that from a Bayesian inference study, 
the results are always finite. 


based on the FIM estimate in the left hand cell, and the 
95% BCI from a Bayesian inference estimate on the right 
hand side. With the FIM results, two things immediately 
stand out : the first is that for sources at relatively low 
redshift (i.e. z < 4), there are some cases where the 
true luminosity distance of the source (represented by the 
dashed black line) is not within the QI. The second is that 
for the majority of sources at z > 6, the error in source 
resolution is > 100%. If we were to stop the analysis here, 
and accept the results of the FIM, we would conclude 
that the eLISA observatory, in its current configuration, 
is only capable of inferring the distance to sources out to 
a maximum redshift of z ~ 6. 

On the right hand side of Fig. 5, the Bayesian inference 
results yield BCIs with finite limits. While the BCI natu¬ 
rally grows in size as we go to higher redshift (due to the 
SNR approaching the threshold for detection), we can see 
that for high 2 sources, we would never claim that they 
were closer than Dl ~ 25 Gpc (or z ~ 2). Furthermore, 
and contrary to the FIM case, the true result is always 
contained within the 95% BCI. We should mention here 
that there is a bias in the retrieved Dl to lower values 
than the true values as we go to higher redshift, which 
we can attribute to the results becoming prior dominated 
at large Dl (low SNR), and our use of prior distribu¬ 
tions that scale as inverse powers of Dl- This bias would 
be removed if we used an astrophysical motivated prior 
distribution with scale parameters determined from the 
ensemble of eLISA detections. 

One final comment on these results is that, in general, 
when the BCIs are largest, this corresponds to a source 


where m\ ~ m 2 . As previously stated, in each of the 
odd waveform harmonic terms, and in some parts of the 
higher even terms, we have factors of 5m = mi — m 2 
appearing. This means that for the almost equal mass 
scenario (which is the case for many of the high redshift 
sources), a lot of the higher harmonic contributions are 
suppressed. This results in a general deterioration in the 
precision of source resolution and an increased error in 
the parameter estimation. 


C. Sky Position 

The second major result in this work concerns the er¬ 
ror in the location of the source in the sky. Once again, 
previous studies based on a FIM analysis have presented 
results where the error box on the sky, A fl is equivalent 
to, or greater than the size of the sky. While GW detec¬ 
tors/observatories are not especially good when it comes 
to sky localisation, one should again expect a finite error 
box for the sky position. 

In Fig. 6 we plot the sky error box for the 120 sources 
using again a FIM analysis (left hand cell) and a Bayesian 
inference analysis (right hand cell). In both cases, the up¬ 
per limit of the plot corresponds to the full area of the 
sky. While not immediately obvious from the plot, the 
results from both methods rarely coincide. Furthermore, 
there is no clear pattern in the results that suggest that 
one method systematically predicts larger uncertainty re¬ 
gions. If we focus on the FIM error box, we can see that 
at z > 9 we start to observe sources where the error box 
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contains the entire sky, and in fact after z > 11 nearly 
all sources in the sample have FIM error boxes that are 
greater than the totality of the sky. On the other hand, 
the Bayesian results are always finite. While the errors 
grow quite quickly as a function of redslrift, the maxi¬ 
mum error box seems to converge to 1 /4 of the total sky 
for distant sources. 

We should once again remind the reader that these 
results reflect the nature of the sources in our sample 
and should not be taken as an absolute reflection of the 
capabilities of the eLISA detector. While we do not ex¬ 
pect the high redslrift behaviour to change, nature may 
be kind enough to give us sources at lower redshift that 
have much smaller error boxes, and allow a coincident 
analysis with EM telescopes. 

D. System Parameters 

In Fig. 7 we plot the 95% QIs and BCIs as a percent¬ 
age error for the two mass parameters M c (row 1) and 
/.t (row 2), and the intervals centered around the median 
subtracted true values (i.e. A/ rMe — X K ) for t c (row 3) and 
l (row 4), for each of the 120 sources. In all cases, the col¬ 
umn on the left hand side represents the FIM estimate, 
while the Bayesian inference estimate is represented by 
the column on the right. To ease interpretation, the y- 
axis for the time of coalescence is given in seconds (with 
the upper and lower boundaries of the plot corresponding 
to ±1 hour), and the y-axis for inclination is given in ra¬ 
dians. We remind the reader that the Bayesian inference 
results are the “gold standard” results, and should be in¬ 
terpreted as being more exact. Finally, we point out here 
that the median value in all cases lies within the BCIs. 

If we first focus on the results for the two mass parame¬ 
ters : in both cases, the estimations are quite remarkable. 
For the FIM analysis, the percentage error in the chirp 
mass is always less than ±1%, and less than ±10% for the 
reduced mass in the vast majority of cases. It is only as 
we approach redshifts of £ ~ 11 that we see the maximum 
error in the reduced mass of ±22%. The Bayesian infer¬ 
ence results paint a similar picture with the vast majority 
of cases returning a BCI of less than ±1% for the chirp 
mass, and ±10% for the reduced mass. Once again, it is 
only at high redshifts that we attain maximum errors of 
±1.2% in M c and ±20% for pt. While there does not seem 
to be much difference in both methods, we should point 
out that in the FIM analysis, most of the the sources at 
redshifts of 2 > 6 have a 100% error in Dl and/or Af2. 
In fact, by referring to these parameters, we should not 
trust any of the FIM results at redshifts of 2 > 6. We 
should once again remind the reader, in order to demon¬ 
strate just how good these results are, that the best ob¬ 
served candidate systems for SMBHBs are at low redshift 
(z < 0.4), yet still have order of magnitude errors in the 


estimation of the mass parameters. 

In row 3, we compare the estimates for the precision in 
the prediction of the time of coalescence. A remarkable 
result of this study is that the maximum Bayesian cred¬ 
ible interval is ±1 hour at z ~ 9. For the vast majority 
of sources, the predictive power of when these systems 
will coalesce is on the order of minutes. Finally, in row 
4 we present the BCIs for the inclination of the source. 
As inclination is important in the estimation of the indi¬ 
vidual masses for an electromagnetically observed binary, 
we can see that it should be possible to make reasonable 
estimations of the inclination of the source out to red¬ 
shifts of 2 ~ 5, with a maximum error in the prediction 
of the inclination of around 60 degrees. 


VI. CONCLUSION 

We have compared the error predictions from both a 
Fisher information matrix and a Bayesian inference ap¬ 
proach for 120 SMBHB sources at redshifts of between 
0.1 <2 <13.2. As has been observed in previous pa¬ 
rameter estimation studies, the FIM once again predicts 
that for certain sources, the error in the estimation of 
luminosity distance and the position of the source in the 
sky is 100%. Taking the FIM results on face value, one 
would conclude that the eLISA detector in its current 
configuration is only capable of source resolution out to 
a redshift of z ~ 6. However, a full Bayesian inference 
analysis shows this conclusion to be false. The error pre¬ 
dictions for the luminosity distance are not only always 
finite, but for each individual source set a minimum dis¬ 
tance below which we would never mistake the source as 
being. It is well known that GW observatories have poor 
angular resolution, and while the error box on the sky 
grows as a function of redshift, this error box is always a 
finite quantity. While not the main goal of this work, we 
also demonstrate that out to redshifts of 13.2, we should 
be able to predict the chirp mass with a maximum er¬ 
ror of < 1%, the reduced mass with a maximum error of 
< 20%, the time to coalescence with a maximum error of 
2 hours, and to a redshift of z ~ 5, the inclination of the 
source with a maximum error of ~ 60 degrees. 

Our study also demonstrated that the posterior dis¬ 
tributions for the parameters describing SMBHBs are 
highly non-Gaussian. This should make one highly sus¬ 
picious of any FIM result as one of the principal condi¬ 
tion hypotheses for usage of the FIM is clearly violated. 
The FIM will continue to be a useful tool when one is 
searching for a quick “order of magnitude” estimate of 
parameter errors. However, we hope that this study has 
clearly shown that for full parameter estimation studies, 
and especially for mission configuration science impact 
studies, other more robust methods are needed. 
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FIG. 7: 95% quantile (left hand side) and credible (right hand side) intervals for M c (row 1), fi (row 2), t c (row 3) and i (row 
4). The intervals are expressed as a percentage error for the mass parameters, while the boundaries of row 3 correspond to an 
error of ±1 hour in the time of coalescence. As we know that the FIM estimates of errors in Dl and AO can be greater than 
100 % at redshifts of z > 6, we should not believe any of the results in the left hand column beyond this redshift as they are 
clearly incorrect also. 
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